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HIGHLIGHTS 


•  A  mathematical  model  of  a  fuel  cell-like  flow  battery  has  been  developed. 

•  The  model  predicts  a  small  flow  penetration  into  a  carbon  paper  electrode. 

•  Reactant  penetrating  an  electrode  limits  the  current  density  of  the  flow  battery. 

•  Predicted  maximum  current  densities  match  published  experimental  results. 
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Flow  batteries  show  promise  for  very  large-scale  stationary  energy  storage  such  as  needed  for  the  grid 
and  renewable  energy  implementation.  In  recent  years,  researchers  and  developers  of  redox  flow  bat¬ 
teries  (RFBs)  have  found  that  electrode  and  flow  field  designs  of  PEM  fuel  cell  (PEMFC)  technology  can 
increase  the  power  density  and  consequently  push  down  the  cost  of  flow  battery  stacks.  In  this  paper  we 
present  a  macroscopic  model  of  a  typical  PEMFC-like  RFB  electrode-flow  field  design.  The  model  is  a 
layered  system  comprised  of  a  single  passage  of  a  serpentine  flow  channel  and  a  parallel  underlying 
porous  electrode  (or  porous  layer).  The  effects  of  the  inlet  volumetric  flow  rate,  permeability  of  the 
porous  layer,  thickness  of  the  porous  layer  and  thickness  of  the  flow  channel  on  the  flow  penetration  into 
the  porous  layer  are  investigated.  The  maximum  current  density  corresponding  to  stoichiometry  is 
estimated  to  be  377  mA  cm-2  and  724  mA  cm-2,  which  compares  favorably  with  experiments  of 
-400  mA  cm-2  and  -750  mA  cm-2,  for  a  single  layer  and  three  layers  of  the  carbon  fiber  paper, 
respectively. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Redox  flow  batteries  have  gained  unprecedented  attention  for 
energy  storage  applications.  Well-suited  for  and  scalable  to 
utility  scale  storages,  they  can  be  employed  to  cope  with  inter¬ 
mittent  power  generation  from  renewable  energy  sources,  such 
as  solar  and  wind  energy  and  are  a  useful  tool  for  load  leveling, 
peak  shaving  and  emergency  back  up  to  improve  the  stability  of 
the  power  grid  [1,2  .  Several  research  groups  are  currently 
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working  on  redox  flow  batteries.  Skyllas-Kazacos  et  al.  [3,4] 
developed  an  all-vanadium  redox  flow  battery  (all-VRFB)  that 
was  motivated  by  earlier  studies  on  iron-chrome  as  a  redox 
couple  [5].  An  advantage  of  the  all-VRFB  is  that  species  crossover 
leads  to  only  a  temporary  capacity  loss  in  contrast  to  a  perma¬ 
nent  capacity  loss  for  the  flow  cells  with  two  or  more  different 
elements  acting  as  anions  or  cations  in  the  electrolyte;  examples 
being  iron-chromium  [6],  bromine-polysulfide  [7]  and  zinc — 
bromine  [8  ,  etc.  Most  recent  studies  are  focused  on  vanadium 
flow  cells,  especially  for  the  all-VRFB.  Savinell  et  al.  [9]  developed 
an  all-iron  flow  cell  with  a  potential  advantage  of  very  low  costs. 
Chiang  et  al.  [10]  invented  a  semi-solid  lithium  rechargeable  flow 
battery  within  using  active  slurry  porous  electrode  material 
through  reaction  chambers  and  an  improved  current  density  is 
achieved.  Aziz  et  al.  11]  has  introduced  an  organic  &  inorganic 
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RFB  with  quinone  reactants  in  an  aqueous  electrolyte.  High  cur¬ 
rent  density  and  discharge  capacity  retention  have  been 
demonstrated. 

Computational  modeling  approaches  have  been  reported  to 
explore  performance  improvements  in  RFB  flow  cells  12-18]. 
Walsh  et  al.  developed  several  mathematical  flow  battery 
models  [14-17]  that  describe  isothermal  14],  non-isothermal 
[15],  hydrogen  [16]  and  oxygen  evolution  [17]  (with  bubble 
formation)  transport  processes  in  electrochemical  reactions. 
These  models  include  diffusion,  convection  and  migration.  A 
dilute-solution  approximation  method  was  incorporated.  Vyn- 
nycky  [18  proposed  an  asymptotic  model  for  a  large-scale  cell 
stack  to  save  computational  time  while  improving  the  calcula¬ 
tion  efficiency. 

Few  experimental  and  numerical  studies  on  details  of  flow 
fields  for  RFBs  have  been  reported.  Kjeang  et  al.  [19]  designed  a 
microfluidic  vanadium  redox  flow  cell  with  a  high  aspect  ratio  Y- 
shape  micro  flow  channel  over  a  porous  carbon  substrate  to 
obtain  a  peak  power  density  of  70  mW  cm-2,  which  is  small  for 
large-scale  energy  storage  applications.  Recently,  Zawodzinski 
et  al.  [20-22]  and  Mench  et  al.  [21,22]  designed  flow  cell  stack 
configurations  for  RFBs  to  improve  current  density  and  power 
density.  These  designs  are  based  on  modern  cell  designs 
employed  in  PEM  fuel  cells  and  allow  operation  at  much  higher 
current  densities  as  compared  to  earlier  conventional  flow 
batteries.  In  these  designs  the  electrolyte  flows  through 
serpentine  channels  behind  a  porous  electrode.  This  is  different 


from  the  conventional  flow  battery  designs  where  the  electrolyte 
flows  through  a  carbon  felt  electrode.  Specifically,  a  polarization 
curve  analysis  method  [20]  was  used  to  analyze  performance 
losses  in  two  types  of  flow  cell  configurations:  (a)  a  simple 
battery  (SB)  with  a  carbon  felt  electrode  embedded  in  the  PVC 
compartments  and  (b)  a  modified  fuel  cell  with  a  serpentine 
flow  channel  fuel  cell  battery  (FCB)  over  the  carbon  felt  (FCB-F) 
or  the  carbon  fiber  paper  electrode  (FCB-P).  A  higher  limiting 
current  density  was  obtained  in  the  FCB-P  than  in  the  SB.  A  peak 
power  density  of  557  mW  cm-2,  which  is  more  than  five  times 
that  obtained  using  conventional  flow  cell  configurations,  was 
achieved  at  60%  state  of  charge  by  adopting  a  “zero-gap” 
serpentine  flow  channel  flow  battery  architecture  [21].  A  supe¬ 
rior  behavior  for  lower  over-potential  was  observed  in  the  flow 
cells  composed  of  two  negative  porous  electrode  layers.  The 
potential  effects  of  thermally  pretreated  porous  electrodes  and 
ion  selective  membranes  on  the  performance  of  the  “zero-gap” 
serpentine  flow  cell  were  also  investigated  [22  .  A  peak  power 
density  of  767  mW  cm-2  was  obtained  by  using  three  layers  of 
carbon  fiber  paper  with  an  oxidative  thermal  pretreatment  and 
the  Nation  212  membrane.  It  was  pointed  out  that  surface 
morphologies  and  oxygen  concentrations  might  alter  the  surface 
functionalization. 

Chen  et  al.  [23]  conducted  experiments  on  a  parallel  flow 
structure  together  with  CFD  analysis  for  flow  distributions  in  the 
all-VRFB,  but  generated  a  very  low  power  density  of 
15.9  mW  cm-2.  It  appears  that  flow  rates  may  be  highly  variable 
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Fig.  1.  The  flow  cell  compartment:  1-end  plates;  2-current  collectors;  3-graphite  plates  with  serpentine  flow  channels;  4-gaskets;  5-porous  electrodes;  6-ion  selective  membrane; 
7-load;  8-electrolyte  tanks  and  9-pumps. 
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over  the  cell,  which  has  motivated  calls  for  an  optimized  flow 
field.  Xu  et  al.  [24]  made  output  performance  comparisons 
among  different  VRFBs  with  no  flow,  serpentine  and  parallel  flow 
channels.  The  VRFBs  with  the  serpentine  flow  channels  was 
found  to  have  a  higher  round  trip  efficiency.  Latha  et  al.  [25] 
carried  out  experiments  on  a  serpentine  flow  field  of  VRFBs. 
The  effects  of  rib  convection,  porous  electrode  compression  on 
flow  distributions  and  pressure  drops  were  studied.  Several 
comprehensive  reviews  of  RFBs  can  be  found  in  the  literature 
[26-29]. 

In  this  paper,  a  model  is  presented  that  helps  explain  the 
maximum  current  density  based  on  rates  of  consumption  of  flow 
reactants  that  can  be  achieved  with  a  single  passage  of  the  “zero- 
gap”  serpentine  flow  channel  flow  battery  design. 

2.  Flow  cell  configurations 

A  redox  flow  battery  cell  with  serpentine  flow  channels  is 
shown  in  Fig.  1  and  this  architecture  is  originally  from  the  work  of 
Zawodzinski  and  Mench  et  al.  [21  .  This  construction  consists  of  the 
following  components:  1-end  plates;  2-current  collectors;  3- 
graphite  plates  with  the  serpentine  flow  channels;  4-gaskets;  5- 
porous  electrodes;  6-ion  selective  membranes;  7-load;  8- 
electrolyte  tanks  and  9-pumps.  The  left  and  right  sides  of  the 
flow  cell  compartment  are  anode  and  cathode,  respectively.  During 
the  discharge  process,  the  electrolytes  in  the  tanks  are  circulated  by 
the  pump  and  pass  through  the  serpentine  flow  channels  over  the 
surface  of  the  porous  electrodes  and  are  accompanied  by  electro¬ 
chemical  reactions. 

3.  Macroscopic  mathematical  model 

A  mathematical  model  is  developed  to  better  understand  the 
mass  transport  mechanism  for  electrolyte  flowing  through  a 
single  passage  of  the  serpentine  flow  channel  and  over  the 
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Fig.  2.  Flow  segment  modeled  of  the  serpentine  flow  channel  with  the  porous  elec¬ 
trode  as  shown  in  Fig.  1. 


porous  electrode.  Particular  attention  is  given  to  the  investigation 
of  the  flow  dynamics  within  the  porous  electrode.  The  transport 
theorem  [30 — 32]  and  averaging  volume  method  [33]  are  applied 
to  establish  the  macroscopic  model,  which  is  largely  based  on 
previous  studies  of  Whitaker  [31-33  ,  Gary  et  al.  [34],  Flowes 
et  al.  [35  ,  Ochoa-Tapia  et  al.  [36,37],  Goyeau  et  al.  [38]  and  Bars 
et  al.  [39]. 

The  flow  passage  segment  modeled  here  is  shown  in  Fig.  2.  The 
current  collector  and  ion  selective  membrane  bounded  by  the 
porous  electrode  and  the  serpentine  flow  channel  are  not  repre¬ 
sented  in  Fig.  2.  The  flow  pattern  could  be  considered  as  a  2D  (X,  Y) 
simplification  of  a  3D  (X,  Y,  Z)  system.  The  2D  macroscopic  model  is 
show  in  Fig.  3.  It  illustrates  macroscopic  physical  structures  for  the 
porous  electrode,  including  a  representative  “average”  elements  for 
the  porous  solids,  adjoining  electrolyte  fluid,  and  the  interface  be¬ 
tween  electrolyte  fluid  and  porous  solids.  Here,  Qf  and  Qp  denote 
the  domain  of  the  flow  channel  and  the  porous  electrode,  respec¬ 
tively.  Xc/,  Zfp  and  Zpm  are  the  corresponding  interfaces  between  the 
current  collector  and  the  flow  channel,  the  flow  channel  and  the 
porous  electrode  and  the  porous  electrode  and  the  ion  selective 
membrane.  For  each  representative  element  for  the  porous  solids 
and  their  neighboring  electrolyte  fluid,  y  and  a  denote,  respectively, 
liquid  phase  and  solid  phase  and  quantities  represent  a  locally 
average  value. 

The  source  of  the  average  dynamic  motion  in  the  porous  layer  is 
derived  based  on  the  assumptions  of  steady,  incompressible, 
Newtonian  and  laminar  flow. 

0  =  -eyV(Py)y  +  flyV2Uy  -  fly£yk^(Uy)  +  SyPyg  (1) 

Eq.  (1)  is  the  macroscopic  mathematical  model  for  the  dynamic 
motion  in  the  porous  layer.  If  the  second  term  (viscous  force 
term)  of  the  right  side  for  Eq.  (1)  is  neglected,  then,  it  reduces  to 
Darcy's  law,  which  states  that  the  superficial  velocity  is  propor¬ 
tional  to  the  permeability  of  the  porous  material  and  the  gradient 
of  the  average  pressure  exerted  on  the  flow  but  inversely  pro¬ 
portional  to  the  dynamic  viscosity  of  fluid.  It's  easy  to  infer  that 
Darcy's  law  deals  with  the  flow  pattern  without  regarding 
viscous  force  term  (Brinkman  term)  [40].  However,  the  viscous 
force  term  is  taken  into  consideration  in  this  macroscopic 
mathematical  model,  which  incorporates  Darcy's  law  and 
Brinkman's  term.  Thus,  Eq.  (1)  is  named  as  Brinkman-Darcy 
model.  Navier-Stokes  motions  are  applied  in  the  flow  channel. 
The  non-dimensional  forms  for  the  dynamic  motions  and 
boundary  conditions  are  shown  below.  In  later  discussions,  ev  fiyy 
py,  ky(T ,  <Py>J  and  <uy>  will  be  replaced  with  s,  jl,  p,  /c,  <Pp> ,  and 
<up>  or  <vp>. 


3.1.  Non-dimensionalization 

Based  on  the  assumptions  mentioned  above  and  neglect  the 
effects  from  gravity,  the  normalized  forms  for  fluid  dynamic 
motions  become  corresponding  Eqs.  (2)  and  (3)  in  X  and  Y  di¬ 
rections  in  the  flow  channel.  Eqs.  (4)  and  (5)  are  expressed  for 
the  dimensionless  flow  patterns  in  the  porous  layer  in  X  and  Y 
directions,  respectively.  Non-dimensional  parameters  are  defined 
in  Table  1. 
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Fig.  3.  2D  schematic  showing  a  macroscopic  diagram  of  flowing  through  the  flow  channel  and  over  the  porous  electrode,  including  a  representative  porous  element:  a  for  solid 
phase,  y  for  liquid  phase,  rpe  denotes  the  average  radius  and  ty  represents  the  characteristic  length  for  y  phase  within  Vya. 
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between  the  flow  channel  and  the  porous  layer  is  the  most  chal¬ 
lenging.  The  earliest  work  on  the  boundary  condition  between  the 
flow  channel  and  porous  medium  was  done  by  Beavers  and  Joseph 
[41  ,  a  semi-empirical  slip  boundary  was  proposed  as  Poiseuille 
fluid  flowing  over  the  porous  permeable  wall.  Further  studies  fol¬ 
lowed  by  Saffman  [42],  Neale  [43],  Williams  [44  ,  Ochoa-Tapia 
[36,37],  Vafai  [45,46],  Alazmi  [47],  Kim  [48],  Poulikakos  [49], 
Chandesris  [50,51]  and  others.  Whether  the  stress  jump  and  con¬ 
tinuity  of  velocity  exist  at  the  interface  are  the  major  questions.  A 
non-physical  effective  viscosity  (neff  =  #l£)  was  proposed.  In  this 
model,  a  common  interfacial  boundary  condition  for  BC#4  is 
adopted  that  accounts  for  the  continuity  of  velocity  and  normal 
stress  [46,48,49,52,53  .  Informative  review  papers  on  BC# 4  are 
found  in  Refs.  [52,53  .  The  BC# 7  as  the  interface  between  the 
porous  layer  and  the  ion  selective  membrane  is  considered  as  “no 
slip”,  which  states  that  the  velocity  at  BC#  7  is  zero  and  it  is  based  on 
the  stationary  ion  selective  membrane. 


Seven  non-dimensional  boundary  conditions  are  described  in 
Fig.  3  and  Table  2.  The  BC#1  is  assumed  to  be  the  ideal  plug  flow 
inlet  rather  the  ideal  parabolic  flow  inlet  and  the  reason  is  that  the 
entrance  electrolyte  flow  is  from  a  perpendicularly  oriented  feed 
pipe  and  this  entrance  flow  is  more  likely  to  be  uniform.  The 
determination  of  BC#4  as  the  interface  boundary  condition 


Table  1 


Table  2 

Non-dimensional  boundary  conditions:  BC#  1,  BC# 2,  BC# 3,  BC# 4,  BC# 5,  BC# 6  and 
BC#  7. 
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Table  3 

Geometry,  electrolyte  properties  and  operation  parameters. 


Parameters 

Value 

Sources 

Thickness  of  flow  channel  (t/) 

0.1  cm 

Measured 

Width  of  flow  channel  (wf) 

0.1  cm 

Measured 

Length  of  flow  channel  (L) 

2  cm 

Measured 

Thickness  of  a  single  porous  layer  ( tp ) 

0.041  cm 

[21] 

Density  of  electrolyte  ( p ) 

1.354  g  cm'3 

[55] 

Dynamic  viscosity  of  electrolyte  (p) 

4.928  x  10'3  Pa  s 

Working  temperature  (T) 

298.15  K 

Initial  ion  concentration  (c) 

0.001  mol  cm'3 

[21] 

Inlet  volumetric  flow  rate  (Qin) 

20  ml  min-1  & 

20  ml  min-1  to 

30  ml  min-1 

[12,21] 

4.  Pre-computational  analysis 

4.2.  Physical  parameters 

Before  employing  a  computational  fluid  dynamics  (CFD)  solver- 
COMSOL  Multiphysics  based  on  discretized  domains,  boundary 
conditions  and  related  algorithms,  2D  physical  parameters  are 
determined.  Parameters  for  geometry,  electrolyte  properties  and 
operating  conditions  are  given  in  Table  3.  The  RFB  flow  cell  with  the 
serpentine  flow  channels  is  evolved  from  direct  methanol  fuel  cell 
(DMFC)  or  proton  exchange  membrane  fuel  cell  (PEMFC)  laboratory 
cell  hardware,  and  the  thickness  (or  depth),  width  and  length  for  a 
single  passage  of  the  serpentine  flow  channel  are  measured  to  be 
0.1  cm,  0.1  cm  and  2  cm,  respectively.  The  thickness  for  a  single 
layer  of  carbon  fiber  paper  is  -0.041  cm  [21].  The  thickness  of  the 
porous  layer  is  altered  from  0.041  cm  to  0.287  cm  (the  number  of 
the  porous  layers  ranges  from  1  to  7)  to  investigate  the  effects  of 
thickness  of  the  porous  layer  on  the  mass  flow  transport  from  the 
flow  channel  into  the  porous  layer.  This  model  is  considered  as 
isothermal  with  an  approximation  as  the  working  temperature  of 
298  I<  54].  The  density  and  dynamic  viscosity  of  the  electrolyte  are 
given  in  Table  3.  Walsh  et  al.  [13-17  ,  Qiu  et  al.  55]  treated  the  bulk 
electrolyte  as  water  while  Knehr  et  al.  [12  found  the  density  and 
dynamic  viscosity  of  negative  and  positive  electrolyte  were  slightly 
different.  The  inlet  volumetric  flow  rate  Qn  =  20  ml  min-1  gives 
U[n  =  33.3  cm  s”1  based  on  the  dimensions  of  the  cross  section  of 
the  flow  channel  and  Qn  is  altered  from  5  ml  min”1 
(ziin  =  8.3  cm  s”1)  to  30  ml  min”1  (uin  =  50  cm  s”1)  to  study  the  inlet 
volumetric  flow  rates  on  the  flow  penetration  into  the  porous  layer. 
Table  4  gives  the  properties  of  the  porous  layer,  including  three 
kinds  of  foametal,  carbon  fiber  paper  and  graphite  felt  porous 
materials.  The  foametal  specimens  were  used  for  the  experimental 
studies  of  Beavers  and  Joseph  [41  ]  to  analyze  the  boundary  condi¬ 
tions  of  the  Poiseuille  flow  over  a  permeable  interface.  The  carbon 
fiber  paper  [28]  and  the  graphite  felt  [24]  are  typical  porous  ma¬ 
terials  for  the  porous  electrodes  used  in  the  RFB  flow  cells. 


Table  4 

Properties  of  porous  layer  materials. 


Porous  layer 
samples 

Pore  (dpo)  or 
fiber (dfi) 
diameter, 
[pm] 

Porosity  e 

Permeability 
k,  [m2] 

Estimated 
permeability 
kck,  [m2]3 

Sources 

Foametal  A 

406  (pore) 

0.78 

9.7  x  10~9 

9.0  x  10'9 

[38,40] 

Foametal  B 

864  (pore) 

0.78 

3.94  x  10  8 

4.1  x  10'8 

Foametal  C 

1140  (pore) 

0.79 

8.2  x  10'8 

8.1  x  10'8 

Carbon  fiber 

10  (fiber) 

0.8 

2.0  x  10'11 

2.3  x  10'10 

[28] 

paper 

Graphite  felt 

17.6  (fiber) 

0.7 

Not  given 

2.1  x  10'10 

[24] 

a  Estimated  by  Eq.  (4). 


Pore  diameters  are  given  for  the  foametal  A,  foametal  B  and 
foametal  C  while  fiber  diameters  are  given  for  the  carbon  fiber 
paper  and  the  graphite  felt.  Goyeau  et  al.  [38]  made  estimations  of 
the  porosity  for  the  foametals  in  Table  4,  and  the  Kozeny-Carman 
model  [56,57]  shown  as  Eq.  (6)  is  used  to  predict  the  perme¬ 
ability  of  the  porous  samples. 


d2e3 

C*c(  1  -  e)2 


where,  d  is  the  pore  diameter  or  fiber  diameter.  The  permeability  is 
directly  correlated  with  pore  diameter  or  fiber  diameter  while 
inversely  correlated  with  the  Kozeny-Carman  constant.  Carman 
pointed  out  that  experimental  C/<c  was  around  180  for  the  pore 
diameter.  Walsh  et  al.  13-17  and  Qiu  et  al.  [55]  considered  C/<c  as 
5.55  for  the  fiber  diameter,  then  the  estimated  permeability  of  the 
carbon  fiber  paper  ( dj\  =  10  pm,  e  =  0.8)  is  around  2.3  x  10”10  m2 
and  this  is  consistent  with  Tomadakis  et  al.  [58]  who  thought  the 
permeability  was  -1  x  10”10  m2  as  a  3D  random  fiber  model  was 
built  for  the  carbon  fiber  paper  with  10  pm  fiber  diameter  and  0.8 
porosity.  However,  Weber  et  al.  [28]  considered  permeability  as 
2  x  10”11  m2  for  the  carbon  fiber  paper  reported  in  fable  4.  In  this 
study,  the  value  of  2.3  x  10”10  m2  is  employed  for  the  estimated 
permeability  for  the  carbon  fiber  paper  with  10  pm  fiber  diameter 
and  0.8  porosity.  The  permeability  of  the  graphite  felt  is  predicted 
from  Eq.  (6)  as  2.1  x  10”10  m2.  The  relationship  between 


Fiber  diameter  (pm) 


Fig.  4.  The  relationship  between  pore  diameter  or  fiber  diameter  and  permeability 
under  specific  porosities  mentioned  in  Table  4  is  curved  through  using  the  Kozeny- 
Carman  model:  (a)  the  foametal:  #-e  =  0.78,  k  =  9.7  x  10  9  m2;  M-e  =  0.78, 
k  =  3.94  x  10-8  m2;  <-e  =  0.79,  k  =  8.2  x  10~8  m2  and  (b)  the  carbon  fiber  paper: 
A-e  =  0.8,  k  =  2.3  x  10~10  m2;  the  graphite  felt:V-e  =  0.7,  k  =  2.1  x  10'10  m2. 
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permeability  and  pore  diameter  or  fiber  diameter  under  a  specific 
porosity  is  shown  in  Fig.  4,  the  permeability  is  predicted  through 
the  Kozeny-Carman  model  with  considering  Qc  =  180  for  the  pore 
diameter  and  Qc  =  5.55  for  the  fiber  diameter. 

Pore  and  fiber  diameters  have  a  great  effect  on  the  permeability. 
Increasing  the  pore  diameter  from  100  pm  to  1000  pm  or  fiber 
diameter  from  10  pm  to  20  pm,  the  permeability  is  elevated  by  a 
factor  of  100  and  4  for  the  pore  diameter  and  fiber  diameter, 
respectively.  Active  surface  area  (as)  of  porous  fibers  can  be  esti¬ 
mated  by  the  filament  analog  model  developed  by  Carta  et  al.  [59] 
within  building  the  porous  fibers  as  cylinder  models  as  shown  in 
Eq.  (7) 


The  variations  of  the  active  surface  area  with  fiber  diameters  are 
shown  in  Fig.  5  for  the  carbon  fiber  paper  and  the  graphite  felt.  The 
fiber  diameter  strongly  affects  the  active  surface  area.  Increasing 
the  fiber  diameter  from  5  pm  to  15  pm,  the  active  surface  area  is 
decreased  by  a  factor  of  3.  However,  the  porous  fibers  in  reality  are 
not  perfect  cylinders,  and  the  complex  shapes  of  the  porous  solids 
are  beyond  our  discussions.  The  tradeoff  between  the  permeability 
and  the  active  surface  area  needs  further  study. 

The  values  for  the  inlet  volumetric  flow  rates,  inlet  velocities 
and  Reynolds  numbers  are  given  in  Table  5.  Because  of  the  Reynolds 
number  ranging  from  22.8  to  137,  the  flow  regime  is  considered  as 
laminar  flow. 

5.  Results  and  discussion 

5.2.  Normalized  flow  patterns 

The  flow  hydrodynamics  are  dominated  by  Navier-Stokes  mo¬ 
tions  in  the  flow  channel  Qf  and  Brinkman-Darcy  model  in  the 


Table  5 

Value  for  Qin,  ifin  and  Re  using  geometry  dimensions  given  in  Table  3. 


Qin  (ml  min  !) 

It  in  (cm  s  1) 

Re 

5 

8.3 

22.8 

10 

16.7 

45.9 

15 

25 

68.7 

20 

33.3 

91.5 

25 

41.7 

115 

30 

50 

137 

porous  layer  Qp.  In  order  to  clearly  see  the  flow  patterns,  larger 
permeability  (9.7  x  10-9  m2)  of  the  porous  layer  is  chosen.  The  flow 

velocities  in  the  X  direction  (ujf  and  (up)*)  are  illustrated  below 
while  the  flow  velocities  in  the  Y  direction  (i^and  (vp)*)  are  not 
shown  because  in  most  of  the  flow  channel  and  the  porous  layer, 
the  velocities  vj  and  (vp)*  are  much  smaller  compared  with  corre¬ 
sponding  and  (Up)*  except  at  the  very  beginning  and  end  of  the 
flow  channel  and  the  porous  layer.  Fig.  6  describes  the  trend  of 
normalized  and  (up)*  along  X*  ranging  0  to  1  for  six  different  Y* 

under  the  ideal  plug  flow  inlet  condition  as  the  Re  (Reynolds 
number)  is  91.5.  In  the  ^/domain,  ujt  climbs  from  1  to  1.4  as  the  flow 

regime  changes  from  uniform  to  fully  developed  state.  However,  at 
near  to  the  front  edge  of  the  Qp  region,  ( up )*  sharply  increases  to  a 
peak  value,  then  decreases  as  the  fully  developed  region  is 
approached.  The  hydrodynamic  entrance  length  Le  is  calculated  to 
be  0.46  cm  as  the  centerline  velocity  is  99%  of  the  fully  developed 
centerline  velocity.  About  75%  of  the  flow  channel  is  fully  devel¬ 
oped.  For  the  case  of  conventional  pipes  or  circular  ducts,  the 
normalized  entrance  length  L*  is  ~6%  of  Reynolds  number  for  the 
laminar  flow.  Atkinson  et  al.  [60]  and  Chen  et  al.  [61]  proposed 
correlations  (Eqs.  (8)  and  (9))  between  the  non-dimensional  hy¬ 
drodynamic  flow  entrance  length  L*  and  the  Reynolds  number  for 
the  laminar  flow  through  the  parallel  plate  flow  channel. 


Fiber  diameter  (pm) 


Fig.  5.  The  variations  of  active  surface  area  (as)  for  the  carbon  fiber  paper  and  the 
graphite  felt  as  the  fiber  diameter  increases  from  1  pm  to  20  pm  by  using  the  filament 
analog  model  [59]  within  building  porous  fibers  as  cylinder  models. 


I*  =  C,  +  C2Re  (8) 

L>(QWTT)  +  C^  (9) 

where,  L*  =  Le/t/,  C\  =  0.625,  C2=C5=0.044,  C3  =  0.63,  C4  =  0.035. 

The  corresponding  estimated  Le  are  0.465  cm  and  0.418  cm  from 
correlations  (8)  and  (9).  These  entrance  length  predictions  have  a 
reasonable  agreement  with  the  numerically  computed  entrance 
length  of  0.46  cm.  Some  back  flow  occurs  in  the  porous  layer  with  a 
small  X*  as  shown  in  Fig.  7.  Using  a  smaller  inlet  volumetric  flow 
rate  (Qm  =  5  ml  min-1)  of  8.3  cm  s-1  for  1%  with  22.8  for  Re,  Le  is 
computed  to  be  0.134  cm.  Thus,  the  extent  of  the  fully  developed 
region  is  over  93%  of  the  flow  channel  (see  Fig.  8).  The  estimated  Le 
from  Eqs.  (8)  and  (9)  are  seen  to  be  reasonable  with  the  numerical 
result  of  0.134  cm  at  0.163  cm  and  0.135  cm,  respectively.  If  Qm  is  in 
the  range  of  20  ml  min-1  to  30  ml  min-1  (Knehr  et  al.  [  12  ),  then  the 
corresponding  values  of  Le  then  lie  in  the  range  of  0.46  cm-0.67  cm. 
This  illustrates  that  roughly  65%-75%  of  the  flow  channel  is  fully 
developed.  In  later  sections,  the  volume  flow  rate  of 
Qin  =  20  ml  min-1  (1%  =  33.3  cm  s-1)  is  used  to  analyze  the 
maximum  current  density  for  the  carbon  fiber  paper  electrodes.  It 
is  emphasized  that  our  model  is  not  the  same  as  a  parallel  flow 
plate  channel.  For  latter,  the  BC#4  should  be  “no  slip”  instead  of 
continuity  of  the  velocity  and  the  normal  stress  at  the  interface 
mentioned  in  Table  2.  Thus,  some  accuracy  is  lost  for  the  entrance 
length  prediction  when  correlations  (8)  and  (9)  are  used.  From  both 
Figs.  6  and  8,  an  over-shoot  is  found  in  the  velocity  profile  of  the 
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Fig.  6.  The  trend  of  normalized  ui  and  (up)*  with  X*  from  0  to  1  at  six  V*:  0.390,  0.319, 
0.284,  -0.071,  -0.142  and  -0.213  in  the  £2/ and  the  Qp  (/<  =  9.7  x  10-9  m2,  e  =  0.78, 
Re  =  91.5,  tf  =  0.1  cm,  tp  =  0.041  cm). 


porous  layer  with  a  small  X*.  This  is  expected,  and  is  caused  by  the 
assumption  of  the  ideal  plug  flow  at  the  entrance  interface 
providing  momentum  to  drive  the  fluid  into  the  porous  layer  near 
to  the  left  wall  (BC# 5).  However,  this  over-shoot  disappears  when 
the  ideal  parabolic  flow  inlet  condition  is  adopted  as  mentioned  in 
the  work  of  Ke  [62].  Fig.  9  describes  the  trend  of  normalized  zi|  and 

(zip)*  along  Y*  from  -0.291  to  0.709  for  four  specific  X*:  0,  0.1,  0.2 
and  0.7  representing  the  entrance  through  the  developing,  devel¬ 
oped  and  fully  developed  regions.  It  is  seen  that  the  normalized 
average  velocity  ul  decreases  from  1  to  0.921,  then  increases  to 

Javg 

0.934  while  (up)*avg  exhibits  an  opposite  trend  and  climbs  from  0  to 
0.179,  then  decreases  to  0.150  as  the  fully  developed  region  is 
approached.  The  normalized  averages  uj?  and  ( up)*avg  for  Y*  withX* 
are  described  by  Eqs.  (10)  and  (11),  respectively. 


7 


U 


Jc 


avg  J  y* 


tf  +  h 
tf 


tf+tp 


u}(X*,Y*)dr 


(10) 


(11) 


Eqs.  (10)  and  (11)  give  the  average  dimensionless  velocity  in  a 
specific  X*  along  Y*. 


Fig.  8.  The  trend  of  normalized  u2  and  ( up )*  with  X*  from  0  to  1  at  six  Y*:  0.390,  0.319, 
0.284,  -0.071,  -0.142  and  -0.213  in  the  Qf  and  the  Qp  ( k  =  9.7  x  10-9  m2,  e  =  0.78, 
Re  =  22.8,  tf  =  0.1  cm,  tp  =  0.041  cm). 
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BC#3 

outlet 


BC#6 


Fig.  9.  Non-dimensionalization  flow  configurations  for  u^  and  (up)*  with  Y*  from  entrance  to  developing,  developed  and  fully  developed  region  in  the  £2/and  the  Qp  at  fourX*:  0, 0.1, 
0.2  and  0.7  (k  =  9.7  x  10~9  m2,  e  =  0.78,  Re  =  91.5). 


5.2.  Inlet  volumetric  flow  rate  &  permeability  effects 

The  effects  from  inlet  volumetric  flow  rates  on  flow  distribu¬ 
tions  for  the  fully  developed  region  are  shown  in  Fig.  10.  The  flow 
penetration  into  the  porous  layer  increases  as  the  inlet  flow  rate 


climbs  from  10  ml  min-1  to  30  ml  min-1.  Fig.  11  illustrates  the  ef¬ 
fects  of  the  permeability  on  the  enhancement  of  flow  penetration 
into  the  porous  layer.  This  clearly  shows  that  the  permeability  has  a 
more  dominant  influence  than  the  inlet  volumetric  flow  rate  on  the 
increase  of  flow  penetration  into  the  porous  layer. 


Fig.  10.  Ufa  nd  <up>  profiles  along  the  Y  direction  for  the  fully  developed  region:  inlet 
volumetric  flow  rate  Q,n  climbs  from  10  ml  min  1  to  30  ml  min-1  (u!n  ranges  from 
16.7  cm  s'1  to  50  cm  s-1,  k  =  9.7  x  10-9  m2,  e  =  0.78,  t/=  0.1  cm,  tp  =  0.041  cm). 


Fig.  11.  Uf  and  <up>  profiles  along  the  Y  direction  for  the  fully  developed  region: 
permeability  increases  from  2.3  x  10  10  m2  to  8.2  x  10  8  m2  for  the  foametal  and  the 
carbon  fiber  paper  porous  material  in  Table  4  (t/  =  0.1  cm,  tp  =  0.041  cm, 
Uin  =  33.3  cm  s-1). 
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5.3.  Thickness  of  porous  layer  &  flow  channel  effects 

The  effects  of  the  thickness  of  the  porous  layer  and  the  flow 
channel  on  flow  distributions  along  the  Y  direction  have  been 
studied.  Fig.  12  depicts  the  flow  distributions  for  the  fully  devel¬ 
oped  region  at  three  different  thicknesses  of  the  porous  layer: 
0.041  cm  (single  porous  layer),  0.123  cm  (three  porous  layers)  and 
0.205  cm  (five  porous  layers)  with  the  fixed  0.1  cm  flow  channel 
depth.  More  flow  is  transferred  from  the  flow  channel  into  the 
porous  layer  with  a  larger  thickness  of  the  porous  layer.  The  in¬ 
fluence  of  the  depth  of  the  flow  channel  is  illustrated  in  Fig.  13 
under  the  same  inlet  velocity  condition.  More  flow  goes  into  the 
porous  layer  with  thinner  flow  channels  even  though  this  corre¬ 
sponds  to  a  smaller  inlet  volumetric  flow  rate. 

5.4.  Volumetric  flow  rate  in  the  porous  layer 

The  flow  penetration  into  the  porous  layer  can  be  estimated  by 
the  following  equation: 

o 

Qp  =  wp  J  «up»fddY  (12) 

-b 

The  relationship  between  Qp  and  tp  under  inlet  volumetric  flow 
rate  Qjn  ranging  from  5  ml  min-1  to  30  ml  min'1  for  the  typical 
carbon  fiber  paper  used  in  the  RFB  flow  cells  is  shown  in  Fig.  14.  The 
value  of  Qp  is  a  small  fraction  of  the  inlet  volumetric  flow  rate.  For 
example,  for  tp  of  0.123  cm  (three  porous  layers),  Qp  is  estimated  to 
be  3.74  x  10'4  cm3  s'1  when  Qin  is  5  ml  min'1. 


Fig.  12.  Uf  and  <up>  profiles  along  the  Y  direction  for  three  different  thicknesses  of  the 
porous  layer:  0.041  cm  (single  porous  layer),  0.123  cm  (three  porous  layers)  and 
0.205  cm  (five  porous  layers)  for  the  fully  developed  region  ( tf  =  0.1  cm, 
Uin  =  33.3  cm  s-1,  k  =  9.7  x  10-9  m2,  e  =  0.78). 


5.5.  Limiting  and  maximum  current  density 

As  electrolyte  flows  through  the  porous  electrode,  the  elec¬ 
trochemical  reactions  at  the  electrolyte  interface  within  the 
porous  layer  result  in  current  flow,  which  is  determined  by  the 
species  flux  to  the  solid  surface,  Faraday  constant  and  the  number 
of  transferring  electrons  per  mole  of  reactant  consumed  during 
the  reaction.  The  species  flux  to  the  porous  solid  interface  is 
related  to  diffusion,  convection  and  migration  processes.  The 
current  density  is  limited  to  the  situation  where  the  concentration 
of  reactant  at  the  solid  electrode  interface  goes  to  zero.  In  this 
case,  the  concentration  over-potential  contributes  a  large  effi¬ 
ciency  loss.  In  order  to  minimize  a  concentration  over-potential, 
the  operating  current  density  is  usually  chosen  to  be  not  larger 
than  25%  of  the  limiting  current  density.  The  limiting  current 
density  (based  on  approximate  membrane  area)  can  be  estimated 
from  the  Fields  law. 


atpnFDc 

8 


(13) 


where,  n  is  the  number  of  electrons  transferred  per  mole  of  spe¬ 
cies  reacted,  F  is  Faraday  constant,  D  is  diffusivity,  c  is  the  bulk 
electrolyte  concentration,  tp  is  the  thickness  of  the  porous  elec¬ 
trode,  ‘a’  is  the  interfacial  electrode  area  per  unit  of  porous  elec¬ 
trode  volume,  and  <5  is  the  thickness  of  diffusion  layer.  Eq.  (13)  also 
can  be  written  as 


ifim  =  atPnFkmc 


(14) 


Fig.  13.  Uf  and  <up>  profiles  along  the  Y  direction  for  three  different  thicknesses  of  the 
flow  channel:  0.04  cm,  0.07  cm  and  0.1  cm  for  the  fully  developed  region 
(tp  =  0.041  cm,  Uin  =  33.3  cm  s_1,  k  =  9.7  x  10~9  m2,  e  =  0.78). 
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Fig.  14.  The  relationship  between  Qp  and  tp  as  Qjn  ranging  from  5  ml  min  1  to 
30  ml  min  1  for  the  typical  carbon  fiber  paper  used  in  the  RFB  flow  cells 
(k  =  2.3  x  10  10  m2). 


where,  km  represents  the  local  mass  transfer  coefficient  at  the 
interface  between  the  porous  electrode  solid  phase  and  the  bulk 
flowing  electrolyte  solution.  Recently,  Zawodzinski  and  Mench 
et  al.  [21]  tested  a  flow  battery  with  the  serpentine  flow  channel 
modified  from  DMFC  or  PEMFC  flow  cell  as  shown  in  Fig.  15(b),  (c), 
and  (d).  Compared  with  the  conventional  RFBs  as  shown  in 
Fig.  15(a),  this  flow  cell  has  a  thin  porous  layer  (typical  -0.04  cm) 


and  it  can  lead  to  a  lower  ohmic  loss  while  enabling  a  higher 
limiting  current  density.  In  the  limiting  case  of  no  electrolyte 
penetration  into  the  porous  electrode,  the  limiting  current  density 
can  be  estimated  by  Eq.  (15)  [63] 

<U/> 

tfDX 


him  =  0.9783 


nFDc 


Using  Eq.  (15)  and  the  following  values  for  inlet  volumetric  flow 
rate  ( Qjn  =  20  ml  min-1),  the  geometry  for  the  single  flow  channel 
(tf  =  Wf  =  0.1  cm,  L  =  2  cm),  bulk  concentration 

(c  =  0.001  mol  cm-3),  diffusivity  (D  =  10-6  cm2  s-1)  and  assuming 
no  electrolyte  penetration,  the  limiting  current  density  is  estimated 
to  be  about  78  mA  cm-2.  However,  this  predicted  value  significantly 
underestimates  the  reported  maximum  current  density 
-400  mA  cm-2  and  -750  mA  cm-2  observed  for  a  single  layer  and 
three  layers  of  the  carbon  fiber  paper  in  Refs.  [21  ,  respectively.  One 
can  conclude  from  these  estimates  that  convection  in  the  porous 
layer  is  important  to  support  such  large  current  densities.  An  es¬ 
timate  of  the  maximum  current  density  can  be  made  by  assuming 
the  reaction  is  limited  by  the  total  flow  reactant  consumption  in  the 
porous  layer,  or  the  stoichiometric  consumption,  and  ohmic  losses 
are  assumed  to  be  neglected.  Applying  a  mass  balance  and  Fara¬ 
day's  law  gives 


Imax 


nFcQp 

WpL 


(16) 


where,  Qp  is  obtained  in  previous  sections  of  this  paper  and  as 
shown  in  Fig.  14  and  wp  is  equal  to  w/.  In  this  estimate,  we  assume  a 
very  large  interfacial  mass  transfer  coefficient  in  Eq.  (14).  The 
estimation  of  maximum  current  density  for  a  number  of  layers  for 
the  carbon  fiber  paper  is  given  in  Fig.  16.  It  can  be  seen  that  as  the 
number  of  the  porous  layer  increases  from  1  to  3,  the  maximum 
current  density  increases  from  377  mA  cm-2  to  724  mA  cm-2. 
Results  for  one  layer  and  three  layers  of  the  carbon  fiber  paper 
compare  favorably  to  the  experimental  results  of  Zawodzinski  and 
Mench  et  al.  21  ]  as  shown  in  Fig.  16.  The  inlet  volumetric  flow  rate 
(0in  =  20  ml  min-1)  in  the  flow  channel  causes  flow  in  the  porous 
electrode  Qp  of  7.81  x  10-4  cm3  s-1  (single  porous  layer, 
tp  =  0.041  cm)  and  1.50  x  10-3  cm3  s-1  (three  porous  layers, 
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ion  selective  membrane 


r;.  porous  electrode  (or  porous  layer) 

Fig.  15.  (a)  A  thick  porous  electrode  for  the  conventional  RFBs;  (b)  a  RFB  flow  cell  with  the  serpentine  flow  channel  modified  from  DMFC  or  PEMFC;  (c)  a  single  passage  of  the 
serpentine  flow  channel  with  a  thin  porous  layer,  cross  section  view  of  (b);  (d)  adjoin  flow  channels  with  a  thin  porous  layer,  cross  section  view  of  (b). 
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Experimental  results,  Zawodzinski  and  Mench  et  al.,  ref.  [21] 
This  numerical  model,  ideal  plug  flow  inlet 
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Number  of  porous  layers  for  carbon  fiber  paper 

Fig.  16.  The  relationship  between  the  maximum  current  density  possible  from  the 
flow  reactant  supply  and  the  number  of  porous  layers  for  the  carbon  fiber  paper 
(Qin  =  20  ml  min  \  c  =  0.001  mol  cm-3, 1  =  2  cm,  k  =  2.3  x  10-10  m2). 


volumetric  flow  rate  in  the  porous  layer  increases  as  the  thick¬ 
ness  of  porous  layer  increases  from  0.041  cm  to  0.287  cm  (the 
number  of  porous  layers  ranges  from  1  to  7)  for  the  carbon  fiber 
paper  electrodes. 

Moreover,  a  model  for  describing  the  maximum  current  density 
is  developed  which  estimates  the  stoichiometric  availability  of 
reactant  in  the  porous  layer  for  reaction.  Under  the  inlet  condition 
(Qin  =  20  ml  min-1  or  uin  =  33.3  cm  s-1),  the  volumetric  flow  in  the 
electrode  layer,  Qp  increases  from  7.81  x  10-4  cm3  s-1  to 
1.5  x  10-3  cm3  s-1  when  the  electrode  thickness  increases  from  one 
to  three  layers  of  carbon  paper.  Under  these  conditions,  the 
maximum  current  density  increases  from  377  mA  cm-2  to 
724  mA  cm-2  as  the  electrode  thickness  is  increased.  Predicted 
maximum  current  densities  are  compared  with  experimental  re¬ 
sults  with  agreement  well  within  the  uncertainty  in  the  measure¬ 
ments  of  limiting  current  densities  and  estimated  values  of  the 
permeability  of  the  porous  layers  (377  mA  cm-2  predicted  vs. 
-400  mA  cm-2  measured  for  1  layer,  and  724  mA  cm-2  predicted  vs. 
-750  mA  cm-2  measured  for  3  layers). 
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tp  =  0.123  cm)  from  Fig.  16,  then  the  maximum  current  density  is 
estimated  to  be  377  mA  cm-2  and  724  mA  cm-2,  respectively. 

The  estimates  of  maximum  current  density  above  are  based  on 
the  assumptions  for  electrolyte  properties,  cell  design  parameters 
and  estimated  permeability  of  the  porous  electrode.  The  estimated 
maximum  current  density  is  very  dependent  on  the  permeability  of 
the  porous  layer.  These  calculations  use  a  permeability  of 
2.3  x  10-10  m2,  which  may  be  large  according  to  Weber  et  al.  [28] 
(2  x  10  11  m2).  Using  a  smaller  permeability  predicts  a  much 
smaller  amount  of  flow  penetrating  into  the  porous  layer,  and  thus 
a  lower  maximum  current  density.  The  pressure  drop  along  one 
flow  channel  and  then  along  the  adjacent  returning  flow  channel 
will  create  a  pressure  gradient  driving  force  that  drives  fluid  into 
the  adjacent  flow  channel  through  the  porous  layer.  In  this  case,  the 
distance  between  the  adjacent  flow  channels  are  short  (see 
Fig.  15(d)),  and  this  by-pass  of  electrolyte  through  the  porous 
electrode  will  supply  reactant  to  the  electrode  similar  to  our  model, 
except  modified  by  a  factor  of  about  10  (flow  channel  length 
divided  by  the  average  distance  between  the  two  adjacent  flow 
channels).  Consequently,  the  estimate  of  the  maximum  current 
density  due  to  the  bypass  flow  again  gives  a  reasonable  estimate  of 
the  maximum  current  density  as  compared  to  experimentally  re¬ 
ported  value  for  the  lower  estimate  of  permeability  of  2  x  10-11  m2. 

Our  model  as  shown  in  Eq.  (16)  is  based  on  the  contact  area 
between  the  flow  channel  and  the  porous  electrode,  and  this  area  is 
smaller  than  the  membrane  area.  But,  the  actual  penetration  of  flow 
into  the  porous  layer  is  likely  under-estimated  since  flow  through 
the  porous  layer  can  extend  beyond  the  width  of  the  channel. 
Furthermore,  the  interdigitated  flow  channels  as  mentioned  by 
Aziz  et  al.  11]  would  direct  all  flow  through  the  porous  layer,  and 
consequently  a  higher  maximum  current  density  could  be 
achieved. 


6.  Conclusions 

A  2D  macroscopic  mathematical  model  is  established  to  cap¬ 
ture  the  dynamic  flow  patterns  in  RFBs  represented  by  a  single 
passage  of  a  serpentine  flow  channel  over  a  porous  electrode.  The 
non-dimensional  average  flow  velocity  in  the  flow  channel  de¬ 
creases  from  the  entrance  to  the  fully  developed  flow  region 
while  an  opposite  trend  is  found  in  the  porous  layer.  The 
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Nomenclature 

as  active  surface  area  (cm-1) 

BC  boundary  condition 

c  concentration  (mol  cm'3) 

d  diameter  (pm) 

D  diffusivity  (cm2  s-1) 

F  Faraday  constant  (96,485  C  mol-1) 

g  gravity  acceleration  vector  (980  cm  s-2) 

g  gravity  acceleration  (980  cm  s-2) 

i  current  density  (mA  cm-2) 

k  permeability  of  the  porous  electrode  (m2) 

kya  permeability  tensor  for  y  phase  and  a  phase  of  the  porous 
electrode  (m2) 

I<C  Kozeny-Carman 

L  length  (cm) 

P  pressure  (Pa) 

<P>  average  pressure  (Pa) 

<Pt>t  intrinsic  average  pressure  within  y  phase  (Pa) 

Q  volumetric  flow  rate  (ml  min-1  or  cm3  s-1) 

T  temperature  (K) 

t  thickness  (cm) 

X  X  direction 

Y  Y  direction 

Z  Z  direction 

u  velocity  vector  (cm  s-1) 

u  X  direction  velocity  (cm  s-1) 

<u>  average  X  direction  velocity  (cm  s-1 ) 

V  volume  of  the  representative  porous  element  (pm3) 

v  Y  direction  velocity  (cm  s-1) 

<v>  average  Y  direction  velocity  (cm  s-1) 
w  width  (cm) 
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Greek  symbols 
7  liquid  phase 

e  porosity 

1  characteristic  length  in  a  representative  porous  element 

(nm) 

Id  dynamic  viscosity  (Pa  s) 

p  density  of  electrolyte  fluid  (kg  nrT3) 

2  interface 

a  solid  phase 

Q  domain 

Subscripts 

avg  average  value 

cf  between  current  collector  and  flow  channel 

e  entrance 

eff  effective 

f  flow  domain 

fi  fiber 

fp  between  flow  channel  and  porous  layer 

in  inlet 

Urn  limiting 

max  maximum 

p  porous  domain 

pe  porous  element 

pm  between  porous  layer  and  ion  selective  membrane 

po  pore 

Dimensionless  number 
Re  Reynolds  number  ( pumtflid ) 

References 

[1]  M.  Skyllas-Kazacos,  C.  Menictas,  in:  Proceedings  of  19th  International  Tele¬ 
communications  Energy  Conference-INTELEC  1997,  Melbourne,  Australia, 
October  1997,  pp.  19—23. 

[2]  M.  Skyllas-Kazacos,  D.  Kasherman,  D.  Hong,  M.  Kazacos,  J.  Power  Sources  35 
(1991)  399-404. 

[3]  M.  Skyllas-Kazacos,  M.  Rychcik,  R.G.  Robins,  A.G.  Fane,  M.A.  Green, 
J.  Electrochem.  Soc.  133  (1986)  1057-1058. 

[4]  M.  Skyllas-Kazacos,  M.  Rychick,  R.G.  Robins,  US  Patent  4,786,567,  USA,  1988. 

[5]  Redox  Flow  Cell  Development  and  Demonstration  Project,  NASA  TM-97067, 
National  Aeronautics  and  Space  Administration,  Washington,  DC,  1979. 

[6]  L.H.  Thaller,  Department  of  Energy,  Washington,  DC,  DOE/NASA/1002-79/3, 
National  Aeronautics  and  Space  Administration,  Washington,  DC,  NASA  TM- 
79143,  1979. 

[7]  R.J.  Remick,  P.G.P.  Ang,  USA  Patent  4,485,154,  1984. 

[8]  K.J.  Cathro,  K.  Cedzynska,  D.C.  Constable,  J.  Power  Sources  19  (1987)  337—356. 

[9]  L.W.  Hruska,  R.F.  Savinell,  J.  Electrochem.  Soc.  128  (1981)  18—25. 

[10]  M.  Duduta,  B.  Ho,  V.C.  Wood,  P.  Limthongkul,  V.E.  Brunini,  W.C.  Carter,  Y.- 
M.  Chiang,  Adv.  Energy  Mater.  1  (2011)  511—516. 

[11]  B.  Huskinson,  M.P.  Marshak,  C.W.  Suh,  S.  Er,  M.R.  Gerhardt,  C.J.  Galvin,  X.  Chen, 
A.A.  Guzik,  R.G.  Gordon,  M.J.  Aziz,  Nature  505  (2014)  195—198. 

[12]  K.W.  Knehr,  E.  Agar,  C.R.  Dennison,  A.R.  Kalidindi,  E.C.  Kumbur,  J.  Electrochem. 
Soc.  159  (2012)  A1446— A1459. 

[13]  A.A.  Shah,  R.  Tangirala,  R.  Singh,  R.G. A.  Wills,  F.C.  Walsh,  J.  Electrochem.  Soc. 
158  (6)  (2011)  A671-A677. 

[14]  A.A.  Shah,  M.J.  Watt-Smith,  F.C.  Walsh,  Electrochim.  Acta  53  (2008) 
8087-8100. 


[15]  H.  Al-Fetlawi,  A.A.  Shah,  F.C.  Walsh,  Electrochim.  Acta  55  (2009)  78—89. 

[16]  A.A.  Shah,  H.  Al-Fetlawi,  F.C.  Walsh,  Electrochim.  Acta  55  (2010)  1125—1139. 

[17]  H.  Al-Fetlawi,  A.A.  Shah,  F.C.  Walsh,  Electrochim.  Acta  55  (2010)  3192—3205. 

[18]  M.  Vynnycky,  Energy  36  (2011)  2242—2256. 

[19]  E.  Kjeang,  B.T.  Proctor,  A.G.  Brolo,  D.A.  Harrington,  N.  Djilali,  D.  Sinton,  Elec¬ 
trochim.  Acta  52  (2007)  4942—4946. 

[20]  D.  Aaron,  Z.  Tang,  A.B.  Papandrew,  T.A.  Zawodzinski,  J.  Appl.  Electrochem  41 
(2011)  1175-1182. 

[21]  D.S.  Aaron,  Q,  Liu,  Z.  Tang,  G.M.  Grim,  A.B.  Papandrew,  A.  Turhan, 
T.A.  Zawodzinski,  M.M.  Mench,  J.  Power  Sources  206  (2012)  450—453. 

[22]  Q,  Liu,  G.M.  Grim,  A.B.  Papandrew,  A.  Turhan,  T.A.  Zawodzinski,  M.M.  Mench, 
J.  Electrochem.  Soc.  159  (2012)  A1246-A1252. 

[23]  J.  Chen,  B.  Wang,  H.  Lv,  Adv.  Mater.  Res.  236-238  (2011)  604-607. 

[24]  Q,  Xu,  T.  Zhao,  P.  Leung,  Appl.  Energy  105  (2013)  47—56. 

[25]  T.J.  Latha,  S.  Jayanti,  J.  Power  Sources  248  (2014)  140—146. 

[26]  Z.  Yang,  J.  Zhang,  M.C.W.  Kintner-Meyer,  X.  Lu,  D.  Choi,  J.P.  Lemmon,  J.  Liu, 
Chem.  Rev.  Ill  (2011)  3577-3613. 

[27]  M.  Skyllas-Kazacos,  M.H.  Chakrabarti,  S.A.  Hajimolana,  F.S.  Mjalli,  M.  Saleem, 
J.  Electrochem.  Soc.  158  (2011)  R55-R79. 

[28]  A.Z.  Weber,  M.M.  Mench,  J.P.  Meyers,  P.N.  Ross,  J.T.  Gostick,  Q,  Liu,  J.  Appl. 
Electrochem.  41  (2011)  1137—1164. 

[29]  P.  Leung,  X.  Li,  Carlos  Ponce  de  Leon,  L.  .  Berlouis,  C.T.  John  Low,  F.C.  Walsh, 
RSC  Adv.  27  (2012)  10125-10156. 

[30]  S.  Whitaker,  Chem.  Eng.  Sci.  28  (1973)  139—147. 

[31]  S.  Whitaker,  AIChE  J.  13  (1967)  420-427. 

[32]  S.  Whitaker,  Transp.  Porous  Media  1  (1986)  3—25. 

[33]  S.  Whitaker,  The  Method  of  Volume  Averaging,  Kluwer  Academic  Publishers, 
Norwell,  1999. 

[34]  W.G.  Gary,  K.  O'Neill,  Water  Resour.  Res.  12  (1976)  148-154. 

[35]  F.A.  Howes,  S.  Whitaker,  Chem.  Eng.  Sci.  40  (1985)  1387—1392. 

[36]  J.A.  Ochoa-Tapia,  S.  Whitaker,  Int.  J.  Heat  Mass  Transfer  38  (1995)  2635—2646. 

[37]  J.A.  Ochoa-Tapia,  S.  Whitaker,  Int.  J.  Heat  Mass  Transfer  38  (1995) 
2647-2655. 

[38]  B.  Goyeau,  D.  Lhuillier,  D.  Gobin,  M.G.  Velarde,  Int.  J.  Heat  Mass  Transfer  46 
(2003)4071-4081. 

[39]  M.L.  Bars,  M.G.  Worster,  J.  Fluid  Mech.  550  (2006)  149-173. 

[40]  H.C.  Brinkman,  Appl.  Sci.  Res.  Al  (1947)  27—34. 

[41]  G.S.  Beavers,  D.D.  Joseph,  J.  Fluid  Mech.  30  (1)  (1967)  197-207. 

[42]  P.G.  Saffman,  Stud.  Appl.  Math.  1  (1971)  93—101. 

[43]  G.  Neale,  W.  Nader,  Can.  J.  Chem.  Eng.  52  (1974)  475-478. 

[44]  W.O.  Williams,  Quart.  Appl.  Math.  36  (1978)  255-268. 

[45]  K.  Vafai,  S.J.  Kim,  Int.  J.  Heat  Fluid  Flow.  11  (1990)  254-256. 

[46]  K.  Vafai,  R.  Thiyagaraja,  Int.  J.  Heat  Mass  Transfer  30  (1987)  1391—1405. 

[47]  B.  Alazmi,  K.  Vafai,  Int.  J.  Heat  Mass  Transfer  44  (2001)  1735—1749. 

[48]  S.J.  Kim,  C.Y.  Choi,  Int.  J.  Heat  Mass  Transfer  39  (1996)  319—329. 

[49]  D.  Poulikakos,  M.  Kazmierczak,  J.  Heat  Transfer  109  (1987)  653—662. 

[50]  M.  Chandesris,  D.  Jamet,  Int.  J.  Heat  Mass  Transfer  49  (2006)  2137—2150. 

[51]  M.  Chandesris,  D.  Jamet,  Int.  J.  Heat  Mass  Transfer  50  (2007)  3422—3436. 

[52]  M.  Ehrhardt,  An  Introduction  to  Fluid-porous  Interface  Coupling,  http://www. 
math.uni-wuppertal.de,  (accessed  May  2014). 

[53]  V.  Laptev,  Numerical  Solution  of  Coupled  Flow  in  Plain  and  Porous  Media 
(Ph.D.  thesis),  University  of  Kaiserslautern,  2003. 

[54]  A.  Tang,  J.  Bao,  M.  Skyllas-Kazacos,  J.  Power  Sources  248  (2014)  154—162. 

[55]  G.  Qiu,  C.R.  Dennison,  K.W.  Knehr,  E.C.  Kumbur,  Y.  Sun,  J.  Power  Sources  219 
(2012)  223-234. 

[56]  J.  Kozeny,  Stizungsber  Akad.  Wiss.  Wien.  136  (1927)  271—306. 

[57]  P.C.  Carman,  Trans.  Inst.  Chem.  Eng.  Lond.  15  (1937)  150—166. 

[58]  M.M.  Tomadakis,  T.J.  Robertson,  J.  Compos.  Mater.  39  (2005)  163—188. 

[59]  R.  Carta,  S.  Palmas,  A.M.  Polcaro,  G.  Tola,  J.  Appl.  Electrochem  21  (1991) 
793-798. 

[60]  B.  Atkinson,  M.P.  Brocklebank,  C.C.H.  Card,  J.M.  Smith,  AIChE  J.  15  (1969) 
548-553. 

[61]  R.  Chen,  J.  Fluids  Eng.  95  (1973)  153-158. 

[62]  X.  Ke,  CFD  Studies  on  Mass  Transport  in  Redox  Flow  Batteries  (M.Sc.  thesis), 
Case  Western  Reserve  University,  2014. 

[63]  J.S.  Newman,  Electrochemical  Systems,  Prentice-Hall,  Inc.,  Englewood  Cliffs, 
1973. 


